************************************************************************************************************
** "Renewalist Christianity and the Political Saliency of LGBTs: Theory and Evidence from Sub-Saharan Afria" 
**  Journal of Politics (2015)  
**  Guy Grossman, University of Pennsylvania. ggros@sas.upenn.edu  
************************************************************************************************************

//note: Figures 1, 2, 4 in the article have been produced by using R. 
//      Please see the R codes in the replication package to reproduce these figures. 


capture log close
clear all
set more off
version 13

*******************************************************************************************
* Figure 3. Pentecostals' Belief and Conservative Views (Left Panel) 
*******************************************************************************************

use "Pew.dta"
tabulate DENOMrec, gen(DENOM)
lab var DENOM7 "Pentecostal"

** Definining independent variables

*1) education
recode educ (9=0) (3=2),gen(eduN)
tabulate eduN, gen(EduAttain)
*2) age
recode q97rec (99=.), gen(age)

** Generating outcome variables 

*1) pray more than once a week
gen pray=Q40==1
*2) pray daily
gen PrayDaily=Q64==1
replace PrayDaily=1 if Q64==2
*3) western music have hurt morality in our country
gen wmusic=Q30==1
*4) How important is religion in your life
gen RegImportant=Q42==1
*5) certain are you about belief in God
gen GodCertain=Q48==1
*6) Relationship with God
gen GodRelation=Q49==1
*7) The Bible is the word of God
gen GodBible=Q57==1
*8) Witness healing
gen Healing=Q67A==1
*9) Received a direct revelation from God
gen Revelation=Q67C==1
*10) Devil
gen Devil=Q67D==1
*11) conflict between being a devout religious person and living in a modern society
gen ModernSociety=Q84==1

global controls Q96 EduAttain3 age b99.income
global perceptions GodCertain RegImportant GodRelation GodBible Revelation PrayDaily Healing Devil
su $perceptions $controls

foreach var in $perceptions{
tab `var' DENOM7  if RELIGrec==1, col
ttest `var' if RELIGrec==1,by(DENOM7)
}  

preserve
foreach y in $perceptions {
tempfile tf`y'  
qui parmby "areg `y' DENOM7 $controls if RELIGrec==1, abs(country)", label level(95) saving(`"tf`y'"', replace) ids("`y'")
}

drop _all

foreach y in $perceptions {
append using `"tf`y'"'
}

keep if parm=="DENOM7"
replace idstr="Witnessed the devil" if idstr=="Devil"
replace idstr="Bible word of God" if idstr=="GodBible"
replace idstr="Absolutely certain in belief" if idstr=="GodCertain"
replace idstr="Personal relationship with God" if idstr=="GodRelation"
replace idstr="Witnessed healing" if idstr=="Healing"
replace idstr="Pray daily" if idstr=="PrayDaily"
replace idstr="Religion very important in life" if idstr=="RegImportant"
replace idstr="Direct revelation from God" if idstr=="Revelation"

sencode idstr, gen(modelid)

eclplot estimate min95 max95 modelid, horizontal /* 
*/title("") xtitle("Coefficient estimate with 95% confidence interval", size(2.4)) ytitle("") /*
*/ scheme(s2mono) graphregion(fcolor(white)) name (IntenseBelief) /* 
*/ eplottype(scatter) rplottype(rcap) /* 
*/ xlabel(-0.04(0.02)0.16,labs(2.4)) xscale(r(-0.04(0.02)0.16)) ylab(#10, labs(2.4) nogrid) /* 
*/ estopts(sort msize(medsmall) ms(O) mc(black)) ciopts(lcolor(black) lw(thin)) xline(0, lpattern(dash)) 

*******************************************************************************************
* Figure 3. Pentecostals' Belief and Conservative Views (Right Panel) 
*******************************************************************************************

foreach var in Q85A Q85B Q85C Q85D Q85E Q85F Q85G Q85H  Q85I{
	gen `var'W=`var'==2
}

global morality Q85AW Q85BW Q85CW Q85DW Q85FW Q85HW Q85IW
su $morality $controls

preserve
foreach y in $morality {
tempfile tf`y'  
qui parmby "areg `y' DENOM7 $controls if RELIGrec==1, abs(country)", label level(95) saving(`"tf`y'"', replace) ids("`y'")
}

drop _all

foreach y in $morality {
append using `"tf`y'"'
}

keep if parm=="DENOM7"
replace idstr="Divorce" if idstr=="Q85AW"
replace idstr="Prostitution" if idstr=="Q85BW"
replace idstr="Euthanasia" if idstr=="Q85CW"
replace idstr="Suicide" if idstr=="Q85DW"
replace idstr="Extra-martial sex" if idstr=="Q85FW"
replace idstr="Abortion" if idstr=="Q85HW"
replace idstr="Homosexual behavior" if idstr=="Q85IW"
sencode idstr, gen(modelid)

eclplot estimate min95 max95 modelid, horizontal /* 
*/title("") xtitle("Coefficient estimate with 95% confidence interval", size(2.4)) ytitle("") /*
*/ scheme(s2mono) graphregion(fcolor(white)) name (MorallyWrong) /* 
*/ eplottype(scatter) rplottype(rcap) /* 
*/ xlabel(-0.04(0.02)0.08,labs(2.4)) xscale(r(-0.04(0.02)0.08)) ylab(#10, labs(2.4) nogrid) /* 
*/ estopts(sort msize(medsmall) ms(O) mc(black)) ciopts(lcolor(black) lw(thin)) xline(0, lpattern(dash)) 


*******************************************************************************************
* Table 1. Determinants of LGBT Saliency (Annual Data)
*******************************************************************************************

use "LGBT_long.dta"
** Defininig panel data 
xtset id year, y

** Assigning global to key independent vars

* 1a. Key variables - Highest circulation newspaper
global keyVarsA "articles_total c.cpolity2##c.crenwshare_pop clegit i.christian_mjr##c.crenwshare_pop" 
su $keyVarsA
* 1b. Key variables - Mean available newspapers
global keyVarsM "marticles_total c.cpolity2##c.crenwshare_pop clegit i.christian_mjr##c.crenwshare_pop" 
su $keyVarsM
* 2. Define controls 
global controls "enr malecrime anyelec chdi claids clwb_gdppc claidcommit_pc clpopulation celf1 Colony2 Colony3 Colony4 Colony5 Colony6 Colony7 aregion2 aregion3 aregion4" 
su $controls

eststo clear

* 1) First Column: Model 1. MB2 (Highest circulation paper)
eststo annual_nb2a: nbreg lgbt_articles $keyVarsA ownership $controls i.year, disp(mean) vce(cluster id) nolog irr 
estadd margins, dydx(crenwshare_pop cpolity2 clegit christian_mjr enr)

* 2) Second Column: Model 2. ZINB (Highest circulation paper)
eststo annual_zinba: zinb lgbt_articles $keyVarsA ownership $controls i.year,inflate(_cons) vce(cluster id) nolog irr 
estadd margins, dydx(crenwshare_pop cpolity2 clegit christian_mjr enr)

* 3) Third Column: Model 1. NB2  (Mean available newspapers)
eststo annual_nb2m: nbreg mlgbt_articles $keyVarsM $controls i.year, disp(mean) vce(cluster id) nolog irr 
estadd margins, dydx(crenwshare_pop cpolity2 clegit christian_mjr enr)

* 4) Fourth Column: Model 2. ZINB  (Mean available newspapers)
eststo annual_zinbm: zinb mlgbt_articles $keyVarsM $controls i.year,inflate(_cons) vce(cluster id) nolog irr 
estadd margins, dydx(crenwshare_pop cpolity2 clegit christian_mjr enr)

#delimit;
estout annual_nb2a annual_zinba annual_nb2m annual_zinbm,
cells(b(star fmt(3)) se(par fmt(2))) stats(ll aic bic N) eform  
label varlabels(_cons constant) starlevels(* 0.10) legend mlabels("NB2" "ZINB" "NB2" "ZINB");
#delimit cr

*******************************************************************************************
* Table 2. Unconditional Marginal Effects: Key Independent Variables 
*******************************************************************************************

//note: if you do not have parmest package installed already, please type 
ssc install parmest 

* 1) Row 1-4: Source: Highest Circulation Model: NB2  
tempfile tf1 tf2 tf3 tf4 
est restore annual_nb2a
margins, dydx(1.christian_mjr crenwshare_pop cpolity2 clegit) post
parmest, format(estimate stderr min95 max95 p %9.3f) list(parm label estimate min95 max95 p,clean noobs) lab saving(`"`tf1'"',replace)

* 2) Row 5-8: Source: Highest Circulation Model: ZINB
est restore annual_zinba
margins, dydx(1.christian_mjr crenwshare_pop cpolity2 clegit) post
parmest, format(estimate stderr min95 max95 p %9.3f) list(parm label estimate min95 max95 p,clean noobs) lab saving(`"`tf2'"',replace)

* 3) Row 9-12: Source: Mean available newspapers: NB2  
est restore annual_nb2m 
margins, dydx(1.christian_mjr crenwshare_pop cpolity2 clegit) post
parmest, format(estimate stderr min95 max95 p %9.3f) list(parm label estimate min95 max95 p,clean noobs) lab saving(`"`tf3'"',replace)

* 4) Row 13-16: Source: Mean available newspapers: ZINB  
est restore annual_zinbm
margins, dydx(1.christian_mjr crenwshare_pop cpolity2 clegit) post
parmest, format(estimate stderr min95 max95 p %9.3f) list(parm label estimate min95 max95 p,clean noobs) lab saving(`"`tf4'"',replace)

*******************************************************************************************
* Figure 5. Marginal Effects of Renewalist Population Share 
*******************************************************************************************

est restore annual_zinba
su cpolity2 if e(sample)
global minpol = r(min) 
global maxpol = r(max)
global skiprenw = (r(max)-r(min))/r(N)
margins, dydx(crenwshare_pop) at(cpolity2 =($minpol (1) $maxpol)) level(90) 
    
#delimit;
marginsplot, scheme(s1mono) ytitle("Marginal Effects") title("Renewalist Christianity Marginal Effect", size(4)) 
level(90) xlabel(-5(5)5) yscale(r(-0.5(.5)1.5))  ylabel(-0.5(.5)1.5) yline(0, lp(dash) lc(red)) 
xtitle("Polity IV (centered)", size(3)) name(fig5) ;
#delimit cr

	
	
	
	
	

	
	
	
	
	
	
	
	
	
*******************************************************************************************
*************************************** APENDIX *******************************************
*******************************************************************************************

// Note: Most of the variables are already defined by previous regressions above. 
//       To reproduce, please execute the above commands first. 



*******************************************************************************************
* Figure 1: Distribution of the Dependnet Variable (Political Saliency Distribution Table)
*******************************************************************************************

qui glm lgbt_articles $keyVarsM $controls, fam(poi) cl(id) cformat(%8.3f)
nbvargr lgbt_articles if e(sample), scheme(s1mono) title("Political Saliency Distribution") xtitle("Number of LGBT-related Articles per Year")

*******************************************************************************************
* TABLE 4: Descriptive Statistics 
*******************************************************************************************

#delimit;
global vars lgbt_articles mlgbt_articles articles_total marticles_total ownership   
christian_mjr renwshare_pop enr Politcomp legit polity2 malecrime anyelec 
hdi laids lwb_gdppc laidcommit_pc lpopulation elf1 
Colony1 Colony2 Colony3 Colony4 Colony5 Colony6 Colony7 
aregion1 aregion2 aregion3 aregion4;
#delimit cr

estpost sum $vars,  listwise 
esttab, cells("mean(fmt(a3)) sd(fmt(a3)) min(fmt(a3)) max(fmt(a3)) count") noobs label

******************************************************************************************
* TABLE 5: Alternative Measure of Political Competition 
******************************************************************************************

global keyM "marticles_total c.Politcomp##c.crenwshare_pop clegit i.christian_mjr##c.crenwshare_pop" 
su $keyVarsM

global keyA "articles_total c.Politcomp##c.crenwshare_pop clegit i.christian_mjr##c.crenwshare_pop" 
su $keyVarsA

* 1) First Column: Model 1. MB2 (Highest circulation paper)
eststo rob_nb2a: nbreg lgbt_articles $keyA ownership $controls i.year, disp(mean) vce(cluster id) nolog irr 
estadd margins, dydx(crenwshare_pop Politcomp clegit christian_mjr)

* 2) Second Column: Model 2. ZINB (Highest circulation paper)
eststo rob_zinba: zinb lgbt_articles $keyA ownership $controls i.year,inflate(_cons) vce(cluster id) nolog irr 
estadd margins, dydx(crenwshare_pop Politcomp clegit christian_mjr)

* 3) Third Column: Model 1. NB2  (Mean available newspapers)
eststo rob_nb2m: nbreg mlgbt_articles $keyM $controls i.year, disp(mean) vce(cluster id) nolog irr 
estadd margins, dydx(crenwshare_pop Politcomp clegit christian_mjr)

* 4) Fourth Column: Model 2. ZINB  (Mean available newspapers)
eststo rob_zinbm: zinb mlgbt_articles $keyM $controls i.year,inflate(_cons) vce(cluster id) nolog irr 
estadd margins, dydx(crenwshare_pop Politcomp clegit christian_mjr)

#delimit;
estout rob_nb2a rob_zinba rob_nb2m rob_zinbm,
cells(b(star fmt(3)) se(par fmt(2))) stats(ll aic bic N) eform 
label varlabels(_cons constant) starlevels(* 0.10) legend mlabels("NB2" "ZINB" "NB2" "ZINB");
#delimit cr

******************************************************************************************
* Figure 3: Marginal Effects (Left Panel) 
******************************************************************************************

** The marginal effects of Renewalist population share displayed as a function of countries' level of dominant religion  
 
graph drop _all	
est restore rob_zinba
su crenwshare_pop if e(sample)
global minrenw = r(min) 
global maxrenw = r(max)
global skiprenw = (r(max)-r(min))/r(N)
margins, dydx(crenwshare_pop) at(christian_mjr =(0 1)) level(90)  

#delimit;
marginsplot, scheme(s1mono) ytitle("Marginal Effects") title("Renewalist Christianity", size(4)) 
level(90) xscale(r(-.5 1.5)) xlabel( 0 1) yline(0, lp(dash) lc(red)) xtitle("", size(4)) 
name(fig3);
#delimit cr

******************************************************************************************
* Figure 3: Marginal Effects (Right Panel) 
******************************************************************************************
	
** The marginal effects of Renewalist population share displayed as a function of countries' level of political competition

est restore rob_nb2a
su Politcomp if e(sample)
global minpol = r(min) 
global maxpol = r(max)
global skiprenw = (r(max)-r(min))/r(N)
margins, dydx(crenwshare_pop) at(Politcomp =($minpol (0.1) $maxpol)) level(90) 
    
#delimit;
marginsplot, scheme(s1mono) ytitle("Marginal Effects") title("Renewalist Christianity", size(4)) 
level(90)  yline(0, lp(dash) lc(red)) xtitle("Gov Majority Share", size(4)) 
name(fig4);
#delimit cr

******************************************************************************************
* Table 6. Mean 3-years 
******************************************************************************************

use "avesample3.dta", clear

encode countryname, gen(id)
tab id
xtset id avesample3

* 1) First Column: Model 1. MB2 (Highest circulation paper)
eststo mean3_nb2a: nbreg lgbt_articles ownership $keyVarsA $controls, disp(mean) vce(cluster id) nolog irr 
estadd margins, dydx(crenwshare_pop cpolity2 clegit)

* 2) Second Column: Model 2. ZINB (Highest circulation paper)
eststo mean3_zinba: zinb lgbt_articles ownership $keyVarsA $controls,inflate(_cons) vce(cluster id) nolog irr 
estadd margins, dydx(crenwshare_pop cpolity2 clegit)

* 3) Third Column: Model 1. NB2  (Mean available newspapers)
eststo mean3_nb2m: nbreg mlgbt_articles $keyVarsM $controls, disp(mean) vce(cluster id) nolog irr 
estadd margins, dydx(crenwshare_pop cpolity2 clegit)

* 4) Fourth Column: Model 2. ZINB  (Mean available newspapers)
eststo mean3_zinbm: zinb mlgbt_articles $keyVarsM $controls,inflate(_cons) vce(cluster id) nolog irr 
estadd margins, dydx(crenwshare_pop cpolity2 clegit)

#delimit;
estout mean3_nb2a mean3_zinba mean3_nb2m mean3_zinbm , 
cells(b(star fmt(3)) se(par fmt(2))) stats(ll aic bic N) eform 
label varlabels(_cons constant) starlevels(* 0.10) legend mlabels("NB2" "ZINB" "NB2" "ZINB");
#delimit cr

******************************************************************************************
* Table 7. Mean 5-years
******************************************************************************************

use "avesample5.dta", clear

encode countryname, gen(id)
tab id
xtset id avesample5

* 1) First Column: Model 1. MB2 (Highest circulation paper)
eststo mean5_nb2a: nbreg lgbt_articles ownership $keyVarsA $controls, disp(mean) vce(cluster id) nolog irr 
estadd margins, dydx(crenwshare_pop cpolity2 clegit)

* 2) Second Column: Model 2. ZINB (Highest circulation paper)
eststo mean5_zinba: zinb lgbt_articles ownership $keyVarsA $controls,inflate(_cons) vce(cluster id) nolog irr 
estadd margins, dydx(crenwshare_pop cpolity2 clegit)

* 3) Third Column: Model 1. NB2  (Mean available newspapers)
eststo mean5_nb2m: nbreg mlgbt_articles $keyVarsM $controls, disp(mean) vce(cluster id) nolog irr 
estadd margins, dydx(crenwshare_pop cpolity2 clegit)

* 4) Fourth Column: Model 2. ZINB  (Mean available newspapers)
eststo mean5_zinbm: zinb mlgbt_articles $keyVarsM $controls,inflate(_cons) vce(cluster id) nolog irr 
estadd margins, dydx(crenwshare_pop cpolity2 clegit)

#delimit;
estout mean5_nb2a mean5_zinba mean5_nb2m mean5_zinbm,
cells(b(star fmt(3)) se(par fmt(2))) stats(ll aic bic N) eform 
label varlabels(_cons constant) starlevels(* 0.10) legend mlabels("NB2" "ZINB" "NB2" "ZINB");
#delimit cr

********************************************************************************************
* Table 8. Placebo Test: Agriculture and Corruption Saliency 
********************************************************************************************

use "LGBT_long.dta", clear  

foreach y in agriculture corruption{
eststo `y'_zinb:  zinb  `y'_articles  $keyVarsA ownership $controls i.year,inflate(_cons) vce(cluster id) nolog irr 
	margins, dydx(1.christian_mjr crenwshare_pop cpolity2 clegit enr)
eststo `y'_nb1:   nbreg `y'_articles  $keyVarsA ownership $controls i.year, disp(constant) vce(cluster id) nolog irr  
	margins, dydx(1.christian_mjr crenwshare_pop cpolity2 clegit enr)
eststo `y'_nb2:   nbreg `y'_articles  $keyVarsA ownership $controls i.year, disp(mean) vce(cluster id) nolog irr  
margins, dydx(1.christian_mjr crenwshare_pop cpolity2 clegit enr)
}

#delimit;
estout agriculture_zinb agriculture_nb1 agriculture_nb2 
corruption_zinb corruption_nb1 corruption_nb1,
cells(b(star fmt(3)) se(par fmt(2))) stats(ll aic bic N) eform 
label varlabels(_cons constant) starlevels(* 0.10) legend mlabels("ZINB" "NB1" "NB2" "ZINB" "NB1" "NB2");
#delimit cr	

********************************************************************************************
* Figure 4. Article Slant by Narrator Type 
********************************************************************************************
use "LGBT 2011 articles slant.dta", clear

keep leader ArticleSlant
gen Slant1=1 if ArticleSlant==1
gen Slant2=1 if ArticleSlant==2
gen Slant3=1 if ArticleSlant==3

foreach var in Slant1 Slant2 Slant3{
recode `var' (.=0)
}

collapse  Slant1 Slant2 Slant3, by(leader)
drop if leader==.

reshape long Slant, i(leader) j(tone) 
lab define tone 1 "Negative" 2 "Neutral" 3 "Positive", modify
lab value tone tone

graph bar Slant, by(leader) over(tone) ytitle("LGBTs-related Articles") scheme(s1mono) note("")
graph export "figures/SlantContentAnalysis.pdf", as(pdf) replace
